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Abstract 

In seasonally ice-covered seas and along the margins of perennial ice pack, i.e. in regions with 
medium ice concentrations, the ice cover typically consists of separate floes interacting with each 
other by inelastic collisions. In this paper, hitherto unexplored analogies between this type of ice 
cover and two-dimensional granular gases are used to formulate a model of ice dynamics at the floe 
level. The model consists of: (i) momentum equations for floe motion between collisions, formulated 
in the form of a Stokes-flow problem, with floe-size dependent time constant and equilibrium 
velocity, and (ii) hard-disk collision model. The numerical algorithm developed is suitable for 
simulating particle-laden flow of N disk-shaped floes with arbitrary size distribution. The model is 
applied to study clustering phenomena in sea ice with power-law floe-size distribution. In particular, 
the influence of the average ice concentration A on the formation and characteristics of clusters is 
analyzed in detail. The results show the existence of two regimes, at low and high ice concentration, 
differing in terms of the exponents of the cluster-size distribution and of the size of the largest 
cluster. 
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I. INTRODUCTION 



Recent global climate change, amplified in the polar regions of the earth, has led during 
the last decades to a substantial reduction of the ice extent, especially in the Arctic 
This trend is strongly coupled to a decrease of the surface area covered with thick, multi- 



year ice, and a corresponding expansion of thin, seasonal ice-cover type 



3|. Whether the 



thinning of the ice cover will eventually lead to permanently ice-free summers in the Arctic 



is still a subject of debate {4]. Undoubtedly, however, it substantially changes the response 
of the ice cover to the external forcing. Among other things, lower mechanical strength of 
thin ice makes it more susceptible to deformation, breaking and divergence. This kind of 
highly mobile ice cover, until recently occurring along the edges of the perennial ice pack 
(a so-called marginal ice zone, MIZ) and in subpolar, seasonally ice-covered seas, has been 
incomparably less intensively studied than the central Arctic ice pack. Global climate trends 
suggest that expanding our knowledge of the dynamics of this type of ice cover will become 
increasingly important in the future. 

In this paper, a term 'medium-concentration ice zone' (MCIZ) is used rather broadly 
to describe a strongly fragmented ice cover type, consisting of clearly separated floes, with 
ice concentrations A < 1, i.e., in a dynamic regime between a free drift {A <^ 1) and a 
compact ice cover (A ~ 1). Within MCIZ defined in this way, sea ice possesses a number 
of properties characteristic of granular materials - assemblies of discrete, macroscopic solid 
particles, interacting with each other via dissipative forces (e.g., friction and/or inelastic 
collisions). Dissipative nature of interactions is responsible for distinctive properties of 
granular materials, different from properties of fluids, solids and gases. For a review of 
collective properties and pattern formation in granular materials see, e.g., Aranson and 
Tsimring 5|. 

Sea ice in MCIZ has all properties required to regard it as a particular subclass of granular 
materials: a two-dimensional (2D) granular gas. It consists of discrete solid fioes (particles) 
immersed in water/air ('interstitial fluids') and moving on the sea surface due to the action 
of external forces of wind, currents etc. The floes can be regarded as rigid and indeformable. 
They collide inelastically, i.e., with a loss of kinetic energy. At medium ice concentrations 
floe-floe collisions contribute substantially to the internal stress in the ice 6-^. Equations 
for the slowly- varying and random component of the floes' velocity {q] are fully analogous to 




FIG. 1. Fragments of Landsat images showing clustered sea-ice floes off the Antarctic coast (at 
I21.8°W, 73.5°S) on 18. Feb. 2008 (a) and 26. Jan. 2008 (b). Source: Q. 



those for dissipative gases [10|, lUl . To the best of my knowledge, those analogies have never 
been explored in sea-ice modeling beyond the study by Feltham who applied granular 
flow theory to formulate a simplified MCIZ model with a composite coUisional and plastic 
rheology, and showed that narrow zones of rapid ice movement along the MIZ edge (ice jets) 
can be obtained as an analytical solution to this model. It has been long recognized that 
collisional rheology is very different from the plastic one, routinely used in sea-ice models. 
However, the available collisional rheology models js-S], including the one used by Feltham, 
are based on unrealistic assumptions of fioes of equal size or very narrow size distributions, 
uniformly distributed on the sea surface. 

With all above-mentioned analogies between MCIZ and granular-gas dynamics in mind, 
there are a number of features that make sea ice unique among typically studied granular- 
gas problems. The most important of those features are, firstly, very wide, scale-free size 
distributions of ice fioes, and secondly, size-dependent equilibrium velocities of fioes under 
a given external forcing. The fioe-size distributions (FSDs) observed in many ocean regions 



see 



13 



17j for recent works), are of power-law type P{d) ~ (i ^ where d denotes the 



floe diameter and > 1. Most recent studies, based on combined satellite and airborne 
data covering a very wide range of floe sizes, showed the existence of two 'regimes' in FSDs 
for large and small fioes, respectively, although the nature of the transition between those 
regimes has not been established yet jl^-[l7|. In this work, we assume for simplicity that 
Or is a constant, and concentrate on the influence of the heavy FSD tail on the patterns of 
motion and cluster formation in MCIZ. 
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This study draws from the analogies between sea ice in MCIZ and granular gases, with 
a goal to develop a theoretical and numerical model of MCIZ dynamics at the floe level, 
and to explore the distinct behavior of sea ice compared to other, well studied granular- 
gas systems. Starting from well established equations of sea-ice motion, it is demonstrated 
that, under certain assumptions, the momentum equation describing the floe's response to 
wind and current can be formulated as a Stokes-fiow problem, with a fioe-size dependent 
time constant and equilibrium velocity. Combined with a hard-disk (HD) model of fioe-floe 
collisions, the Stokes equations are used to formulate a numerical model of MCIZ dynamics, 
in which, contrary to previous attempts, fioe-fioe interactions are simulated directly. In the 
further parts of this paper, the model is applied to study one of the fundamental aspects 
of granular-gas dynamics, namely formation of clusters due to irreversible loss of kinetic 
energy of colliding inelastic particles (see II Al below). The general idea is that clusters 



of fioes, observed in many real situations (Fig. [H see also [17|), can be explained as an 
inherent property of sea ice, and not necessarily as a 'passive' result of the external forcing 
(convergence of surface currents). 

A number of possible further applications of the model developed here are beyond the 
scope of this paper and will be explored elsewhere. They include formulation of internal- 
stress parametrization schemes suitable for MCIZ (at present, as mentioned above, only very 
crude collision-rheology models are available and in la^e-scale sea-ice models the internal 
stress is taken into account only at A ^ 1, see, e.g., [18]), and improved formulae for the 
atmosphere-ocean momentum and heat transfer in MCIZ. Some possible consequences of 
cluster formation in sea ice for the above processes, as well as the existence of possible 
feedbacks, are discussed at the end of this paper. 

Apart from sea-ice modeling applications, the results of this work may be also of interest to 
researchers studying behavior of granular gases (particle-laden flows in particular) with wide 
)article-size distributions, which are still relatively unexplored theoretically and numerically 



19|. In particular, hitherto studies of clustering processes in polydisperse granular media 



concentrated either on binary or narrow size distributions |20l-|22j . The model and code 
developed here can be easily applied to other granular-gas systems consisting of particles 
with any given size distribution and size- dependent response to forcing. 

The structure of this paper is as follows: further parts of the Introduction provide a brief 
theoretical background of cluster formation in granular gases, measures of cluster properties 



and methods of molecular-dynamics (MD) modeling of granular gases. Readers familiar with 
physics and modeling of granular materials may skip this part. In Section [Tll the equations 
of ice-floe motion are formulated in a way suitable for MD simulation. The model and its 
numerical implementation are described in Section IIIH followed by the presentation of the 
modeling results in Section [IVl The paper concludes with a discussion in Section |Vl 



A. Clustering in granular gases — mechanisms and cluster characteristics 



In granular gases, clusters - regions with higher than average particle density - are ini- 
tiated by random density fluctuations. Increased collision frequency in regions with higher 
density leads to a local decrease of the 'granular temperature' (kinetic energy associated 
with the random component of the particle motion, analogous to the thermodynamic tem- 
perature of gases), and thus to a local decrease in pressure. As a result, clusters tend to 
grow by attracting particles from their neighborhood, so that they may persist for a long 
time |lo|. This phenomenon is known from real dissipative granular systems and has been 



extensively studied numerically (e.g., |2l|-|27j). The tendency to form clusters, as well as 



their characteristic size, depend on the average properties of the system (particle density, 
coefficient of restitution, particle size and mass distribution, etc.), as well as on the external 
forcing. Most studies on granular materials consider equal-size and equal-mass particles. In 
the case of continuous size distributions, larger particles tend to accumulate in the central 
regions of clusters (a 'self-sorting mechanism', 22|). 

A number of parameters can be defined for the purpose of a quantitative description of 
clustering processes. The most straightforward way of defining clustered/diluted regions is 



in terms of the local partic 



analyzed region (e.g. 



22 



e concentration A relative to the average concentration A in the 



27|). In this paper, another, distance-based approach is taken. 



i.e., the particles are assigned to clusters in such a way that ever y tw o particles separated by 



a distance smaller that dm belong to the same cluster (e.g., [25|, |26|). Cluster properties are 
defined as averages and/or distributions of the properties of particles forming them, see Sec- 
tion llVi Additionally, a frequently used global measure of the spatial particle distribution 
is a radial distribution function grdf{x), describing the expected particle concentration at a 



distance x from the center of a randomly selected particle, normalized with A |20l . |22| . |27] . 



This function is particularly suitable for monodisperse systems 



20[ | , as it can provide a mea- 



sure of characteristic cluster size. For polydisperse systems, floe-floe distance distribution 
gffd{x) proposed in Section HV Bl seems more suitable. 



B. Modeling of granular gases 

At the core of an MD simulation is the detection of time and partners of collisions among a 
(usually very large) number of particles constituting the analyzed system. Hence, algorithms 
for solving typical fluid-dynamics problems, in which the time is advanced by prescribed, 
constant time steps, are unsuitable. Instead, event-driven algorithms (EDA) must be used 



28 



29| . The choice of a proper EDA depends on (i) the character of motion of the particles 
between collisions, i.e., their interactions with the surrounding fluid, and (ii) the nature of 
particle-particle collisions. Most relevant for this work are algorithms designed for a subclass 



of dissipative-dynamics problems called particle-laden flows [29|, in which the motion of the 
fluid is given, i.e., it influences, but is not influenced by particles immersed in it. 

Considering the material properties of sea ice it is justifled to use a hard-disk (HD) 
collision model, assuming that at collision the momentum is transferred instantaneously 
along the line joining the centers of colliding particles 28|]. In other words, collisions are 
pairwise and inflnitely short. EDAs for HD models are built based on event lists, in which 
times of collisions between all pairs of (potentially colliding) particles are stored and updated 
every time a collision takes place or particles' velocities change. Available EDAs differ mainly 
in terms of how the event lists are updated, which may have aprofound influence on the 



computational efficiency of those EDAs. For examples, see 30, |3l| and a review by [29 1. 



II. ASSUMPTIONS AND EQUATIONS 

Let us consider a set of disk-shaped, non-overlapping sea-ice floes with diameters di 
(z = 1, . . . , A^), equal heights h and constant density p, so that their mass centers coincide 
with their geometric centers. Due to the action of external forces, the floes move on the sea 
surface along trajectories Xj(t), where Xj G ^ denotes the position vector within an analyzed 
region ^ C and t denotes time. The state of each floe at time t is given by (xj(t), Uj(t)), 
where u, = dxj/dt denotes the floe's translational velocity. Without interactions with 
neighboring floes, each floe moves independently and its acceleration can be determined from 
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its state and the external forces, as described below. If the trajectories of floes intersect, 
they collide inelastically, with a loss of kinetic energy. 



Floe motion between collisions 



Between collisions, the motion of the i-th floe satisfies the momentum conservation equa- 



tion (e.g., |32|): 

m, + /k X u,^ = Ff,,dV + ^ Fs.dS, (1) 

where k = [0, 0, 1], / denotes the Coriolis parameter, rrii = phSi is the floe mass. Si = T^df/A 
- its upper/lower surface area, and F^^j, j denote the sum of body and surface forces, 
respectively, acting on the analyzed fioe. For simplicity, we will further assume that / = 
and Fb^i = (in particular, the force resulting from the gradient of the geopotential height 
of the sea surface is not taken into account). The net surface force F^^j results from a sum 
of four terms (Fig. [2]): stress acting on the upper and lower surface of the fioe (atmospheric 
and oceanic skin drag, Tha,i and Th^^i, respectively), and pressure acting on the floe's edges 



(atmospheric and oceanic body drag, Tva,i and Tvw,i-, respectively), see, e.g., 

laJQ- The 

vertical surface area exposed to Tya,i and T^,^_j equals dihf and diih — hf), respectively, 
where hf = h{p^ — p)/Pw denotes the floe's freeboard and denotes water density. (More 
realistically, the surface areas depend on the height/depth of ridges/keels, respectively, but 
this effect is not taken into account here.) Thus, ([T]) can be rewritten as: 

m.^ = ^^(^^- + ^^-^) + ^/^^^-.^ + - hMr.^.- (2) 

In the following, the four forcing terms are parameterized with simple drag formulae typically 
used in sea-ice and ocean modeling. The atmospheric forcing depends quadratically on the 
wind speed u^: 

'^ha,i PaC*/ja I I Uq, Tva,i PaC'tia | Ua | Ua; (3) 

where pa denotes the air density. The oceanic forcing depends on the floe velocity relative to 



the sur: 



one |32| : 



ace current u^. A quadratic relationship analogous to ([3]) is replaced with a linear 



'^hw,i PwChwiy^w '-Ij)) '^vw,i PwC^vw{,^w ^i) ■ (4) 

The drag coefficients Chai C^a, Chw, and C^^ are assumed constant. Linearization of 
Thw,i and Tyyj^i has important consequences for numerical algorithms used in the model, as 



wind 













u 


h 





current 



FIG. 2. A side view of a disk-shaped ice floe under the action of wind and current. The arrows 
are not to scale. See the text for explanation of the symbols. 



it enables to formulate thegoverning equation ^ in the form mathematically analogous to 



the Stokes-flow problem 



29| . Namely, from 



we have: 



du,; 



(5) 



where the time constant Tj is given by: 



Ti 



PwCi 



hw 



+ 



4a, 



(6) 



ph ndi 

After enough time and without interactions with other floes, each floe would reach its equi- 
librium (duj/dt = 0), free-drift velocity u^g^i, given by: 



where: 



Ci — Ti 



PaCha _^ 4:{pw - p)PaCv 



(7) 



(8) 



ph TTPwpdi 

As can be seen, due to the form-drag effects Xj and u^g^i are floe-size dependent, as shown 
in Fig. [3] for a selected range of model parameters. 



B. Floe— floe collisions 



Two floes, i and j, moving along intersecting trajectories, collide at time tc which is the 
smallest positive root of the equation expressing the contact criterion for those floes: 

||xi(g -x,(g|| = {di + dj)/2. (9) 

As already mentioned, it is assumed that collisions are inelastic and inflnitely short, and 
that they can be resolved via a HD model. If Uj, uj denote velocities just before collision, 

8 



(a) (b) 




10° 10' 10' 10° 10' 10' 

d. (m) d. (m) 



FIG. 3. Time constant r (a), coefficient C and tlie equilibrium floe velocity \ueq\ for wind speed 
|ua| = 10 m/s (b) in function of the floe diameter di, for a set of ice thickness values h. The lines 
are drawn every 0.5 m from h = 0.5 m to /i = 3.0 m. Thick lines correspond to the model setup 
used in the simulations described in section |IV] (see Table H]). 

the post-collision velocities, Uj and Uj, are given by: 

Uj = Uj — nijiiij, Uj = Uj + mjiijj, (10) 

where: 

Uij = ^ ^ (kij ■ (uj - Uj)) kij, (11) 

irLi + nij 

and kjj is a unit vector pointing from Xj to Xj. The restitution coefficient e G [0, 1] is 
assumed constant. 

III. THE MODEL 

Like other MD models (see Section II Bp . the one developed in this study - fancifully 
called the Small Floe Collider (SFC) - is based on an event-driven algorithm, suitable for 
particle-laden flows. SFC is capable of simulating a set of disc-shaped particles ('floes') 
with an arbitrary size distribution, within a square domain of side length L, with periodic 
boundaries. In order to simplify the treatment of periodic boundaries, during initialization 
the model domain is scaled to a region [—1/2,1/2] x [—1/2,1/2]. All spatially-dependent 
variables (dj, |uj| etc.) are scaled accordingly. In order to speed up the simulations, the 
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model domain is divided into Nc square computational cells in order to reduce the search 
region for potential collision partners of a given fioe. Consequently, two event types must 
be handled: floe-fioe collisions (FFC) and virtual- wall collisions (VWC), corresponding to 



a transfer of a fioe to a neighboring cell 



29|. 



Before we describe the SFC algorithm, two issues are worth discussing. Firstly, when 
discretizing the model equations, it should be noted that the fioes' trajectories between 
collisions can be determined exactly from equation ([5]). With initial conditions Uj(to) = Uj^o 
and Xi(to) = Xj^o, we obtain easily, for At = t — to: 

Xi(t) = Xi,o + AtUeg^i - Ti{ui^o - Ueg,i)(e-^*/^' - 1). (12) 



Nevertheless, after inserting the floe positions given by f[T2|) into the collision criterion (|9]), 
the resulting equation cannot be solved explicitly for tc. Hence, ( !T2l) must be replaced with 
an approximate formula, allowing for a consistent calculation of fioes' positions and collision 
times. In SFC, the following linearly implicit integrator is used: 

Xi(t) = Xi,o + AtUi^o, (13) 

Ui{t) = Ui,o + At/ri{Ueq^i - Ui{t)) . (14) 



It can be shown easily that using (1131) amounts to replacing the exponential term in ( 1121) 
with the first two terms in its series expansion (e^ ~ 1 + x). From (1121) and (fT3!) . the error 
of floes' positions estimation can be calculated, and thus the maximum acceptable Atmax at 
which velocities must be updated. The error tends to be larger for smaller fioes (because of 
both smaller Tj and, typically, larger |uj — Ugg^il, see further section HVT) . 

The second issue requiring some comment is related to the typical numerical problem of 
hard-sphere models: the inelastic collapse due to diverging collision rates between densely 
packed, almost touching particles. In SFC, it is accounted for by the method of Luding 
and McNamara 36(]: if a given fioe participates in more than one fioe-fioe collision within 
a certain time period Ate, corresponding to a collision duration, the restitution coefficient e 
for those collisions is set to 1. In the simulations described in this paper. Ate = 10~^ s has 
been established experimentally as an optimal value. 

The main part of the computational algorithm (Fig. S]) can be summarized as follows: 
At every iteration, the type of the nearest event, its partner(s) and time t^ are determined. 
Subsequently, the model time t is advanced by At = tc — t, the fioes are moved to their new 
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find type, time and partners of the next event 

increase f and Af„pjj,e by Af = - f 

move floes forward (eq. 13) 
and apply periodic boundaries 

compute collision dynamics 
(FFC: eq. 10,1 1 ; VWC: update cell info) 




compute new velocities (eq. 14) 



fully reset event lists within cells 



reset event lists only for 
partners of the last event 




FIG. 4. Flowchart of the SFC model. 



positions, and the event is handled, accordingly to its type (FFC or VWC). In the final part 
of the main loop, the event lists are either updated only for the partners of the last event, 
or fully recalculated for all particles - depending on whether particle velocities have been 
updated at the current time step. 



IV. MODELING RESULTS 

In this section, the results of SFC simulations will be presented, with emphasis on the 
infiuence of the ice concentration A on cluster formation. A detailed analysis of model 
behavior in the full range of parameters (|ua|, e, etc.) is beyond the scope of this 
study. The general results in terms of fioes' motion and cluster formation are described in 
Section IIV A[ A more detailed analysis of the cluster properties and their size distribution 
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is given in Sections IIVBI and IIV CI 



The model configuration is summarized in Table HI In the simulations, the densities 
(p, pw, Pa) and air-ice drag coefficients (Cha, C^a) were set to their typical values found 
in literature (e.g., 32|). The linear ice- water drag coefficients were estimated from their 



quadratic version 3.5 • 10"'^) and equations (l4l[6] -[8|) for |ui| ~ |ueg^j|, assuming Cha = C^a 
and Chw = Cvw The surface current speed relative to the wind speed was estimated with 
the Ekman model 3^. For the ice floes properties, h, d and a^-, typical values observed in 



15| were chosen. The mean ice concentration 



the Weddell Sea during the ISPOL experiment 
A was varied between the subsequent model runs in order to investigate its influence on the 
modehng results. 



In all simulations, = [mq, 0] and the number of floes = 3000. (Larger system sizes 
were tested for A = 0.6 and A = 0.8, but no differences in terms of statistical properties of 
the results were recorded.) The floe radii di, i = 1, . . . ,N, were obtained with maximum- 
likelihood estimation for a power-law distribution with a given exponent (see appendix I A II 
for details). In each case, the size of the model domain was adjusted to the generated set 
of the floe diameters, to obtain the desired average ice concentration A. For each set of 
model parameters, N^ns = 5 ensemble model runs were performed, differing in respect of 
initial conditions, i.e., the floes' initial velocities (drawn from a normal distribution with 
zero mean and standard deviation 1 ■ 10~^ m/s) and random positions. The optimal number 
of computational cells Nc, the time step for updating the floes' velocity At^ax, and the 
collision duration Ate were established experimentally as 10^, 1.0 s and 0.1 s, respectively. 
For Atmax and At^, the largest values that had no significant influence on the modeling 
results were chosen. 



In all cases, the model was run until a quasi- stationary state developed, characterized by 
insignificantly small temporal changes of such global variables as the total kinetic energy, 
collision rate and pressure. Comparisons between the ensemble model runs have shown that 
the initial conditions had no influence on the final values of those variables, although they 
may have influence on the details of the time evolution of the system towards the final state. 
The required number of iterations was established experimentally to 5 ■ 10^ FFC events. 
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(a) 




FIG. 5. Snapshots (floe positions and velocities, in m/s, shown in gray scale) in a clustered state 
for A = 0.6 (a), A = 0.725 (b) and A = 0.85 (c). For better clarity, only equal-size fragments of 
the model domains are shown. 13 



TABLE I. Physica,! cLiid numerica.1 


model parameters 


used in the SFC simulations 




Parameter 


Symbol 


Value 


Units 


ice density 


P 


910 


kg/m3 


water density 


Pw 


1025 


kg/m3 


air density 


Pa 


1.23 


kg/m3 


air-ice skin-drag coef. 


Cha 


1.7 • 10-3 




air-ice form-drag coef. 


Cva 


1.7 • 10-3 




water-ice skin-drag coef. 


Chw 


1.3 • 10-3 


m/s 


water-ice form-drag coef. 


Cvw 


1.3 • 10-3 


m/s 


wind speed 


Wa\ 


10 


m/s 


surface current speed 




0.01|uJ 


m/s 


ice thickness 


h 


1.5 


m 


mean ice concentration 


A 


0.6-0.9 




restitution coef. 


e 


0.85 




FSD slope 


ar 


1.5 




mean floe diameter 


d 


2.0 


m 


No. floes 


N 


3000 


— 


No. computational cells 


Nc 


102 




time step for vel. updates 




1.0 


s 


collision duration 


Ate 


10-1 


s 


No. ensemble model runs 




5 





A. General properties of the solution 

Figure |5] shows selected snapshots of the floes' positions and velocities from three model 
runs with different A. Contrary to the random initial state of the system, in the quasi- 
stationary regime the motion and distribution of the floes is patterned, with clearly developed 
clusters. Differently than in freely cooling granular gases, however, in the forced system 
under study the formation of clusters is strongly related to nonuniform equilibrium velocity of 
floes with different sizes. Larger floes tend to catch up with the smaller ones that accumulate 
in front of them and then slide along their sides, producing characteristic wakes. They are 
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visible very well especially at small ice concentrations (Fig. |5^), but present even in densely 
packed ice-floe fields as narrow, elongated areas of open water upwind from large floes 
(Fig.Et). Within clusters, the floes tend to stay in contact with each other. The translational 
velocity (i.e., velocity averaged over a certain number of collisions) is almost constant for all 
floes building a cluster. This translational motion is disturbed only if a cluster collides with 
another floe/group of floes. Unless such collision leads to a break-up of a cluster, velocities 
of floes constituting it remain highly correlated. As a result, the along- wind velocity of the 
largest floes within clusters is typically lower than their equilibrium velocity (Fig.|6^,c). This 
difference is larger for larger A, i.e., higher collision rates. To the contrary, the smallest floes 
tend to move up to 30% faster than their equilibrium velocity, with a random component 
of velocity in the order of 1-3-10"^ m/s, comparable to observed values, e.g. from the 
MIZEX experiment g], 7|. Understandably, the higher the ice concentration, the lower the 
differences between the largest and the smallest floes (Fig. [Gt). Even though the differences 
|uj — Ueg,j| may be very large for individual floes, in all analyzed cases (i.e., independently 
of A) the total momentum ^Uiimi develops in time towards ^Ueq,iTni, which corresponds 
to floe-mass-weighted average velocity of ~0.23 m/s. At high ice concentrations, this value 
roughly corresponds to the highest density of points in the along-wind velocity scatterplot 
(and to the velocity of the largest floes. Fig. 



B. Measures of clustering 

Similarly as in observed clusters of ice floes (Fig. [1]), the behavior of the system is dom- 
inated by the largest floes (Fig. [5]). It is not surprising considering that among = 3000 
power-law distributed floe diameters with = 1.5, the largest floe occupies almost 10% 
of the total floe area Ftot = ALF' = df, and the ten largest floes occupy over 46% 

of Ftot- This influence is clearly seen in the shape of the radial distribution functions of 
floes Qrdf (see Section II A|) . For the initial and final states of model runs with A = 0.6 
and A = 0.9, grdf is shown in Fig. [71 As can be seen, the dominating feature of grdf is a 
pattern of shallow minima separated with narrow spikes corresponding to the radii of the 
largest floes, i.e., reflecting small floes densely packed along their edges. Interestingly, for 
A = 0.9, the number of small floes surrounding a very large one is roughly proportional to 
its perimeter; hence an approximately linear dependence of the height of the peaks of Qrdfix) 
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FIG. 6. Scatterplots of instantaneous floe velocity components parallel (a,c) and perpendicular 
(b,d) to the wind direction, as a function of the floe radius r, for A = 0.6 (a,b) and A = 0.9 (c,d). 
The lines in (a) and (c) mark the equilibrium velocity \ueq\. 



on X in Fig. [Tt. 

The grdf properties described above are present already in the initial, random floe distribu- 
tions (in Fig.[7t the two curves are hardly distinguishable) and result simply from geometrical 
constraints of non-overlapping floes. Hence, Qrdf is a poor indicator of the degree of cluster- 
ing, especially for high ice concentrations (Fig. [Tb). An alternative measure of clustering, 
suitable for very wide particle-size distributions, is therefore desirable. A straightforward 
candidate is a fioe-fioe distance distribution gfjd, where the fioe-fioe distance is measured 
not between the fioes' centers, but between their edg instead 
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(a) (b) 




Distance {m) Floe-floe distance (m) 

FIG. 7. Radial distribution function Qrdf (a,c) and floe-floe distance distribution ^//d (b,d) at the 
initial (random) and 'final' (clustered) phases for A = 0.6 (a,b) and A = 0.9 (c,d). Symbols along 
the upper border of (a) and (c) show radii of the eight largest floes in the sample. 



of x = ||xj — Xjll, as in Qrdf)- For random floe distribution in space, Qffd is approximately 
distance-independent at small ice concentrations (Fig. [7)d) and slightly increases toward zero 
in a densely packed ice field, again because of the above-mentioned geometrical constraints 
(Fig. [T]!). However, in a clustered state gfjd increases rapidly (up to an order of magnitude) 
towards zero. 

In the next section, the value of dm = 0.5 m was used to define clusters (see Section [L^ . 
Other tested values from the range 0.1-1.0 m gave very similar results. 
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C. Cluster size distribution 



A deeper insight into cluster formation and properties can be gained from an analysis of 
cluster-size distribution. Let Nc denote the number of clusters in the analyzed snapshot, and 
Ck - an nfc-element set of floes belonging to the A;-th cluster {k = 1, . . . , and Uk = 

N). In the following, a total surface area of floes building the k-th cluster, Fc^k, is used as a 
measure of its size: 

Fc,k = ^J2dl (15) 

Because of generally very irregular shapes of clusters, especially at medium ice concentrations 
(Fig. [5]), this quantity seems more appropriate than e.g. the area of the sea surface occupied 
by a cluster, which is difficult to estimate reliably (in particular, a straightforward convex- 
hull approach turned out inappropriate). Obviously, -^c,A; = ^tot = AL"^. The effective 
cluster radius is given by Ve^k = {Fc^k/T^Y^"^ ■ 

Figure |8] shows the rank-order distribution (see Appendix lA 21) of in the final, quasi- 
stationary state for three values of A: 0.6, 0.75 and 0.9. The distributions have two domi- 
nating features. Firstly, over most of the range of values they are of power-law type: 

P(re;A) ~r-i-"(^). (16) 

Secondly, the sizes of the largest clusters clearly deviate from the power-law regime, the 
more so the larger the ice concentration. 

The exponents a of probability density functions P{re), estimated with a maximum 
likelihood method (equation IA2I in Appendix IA2p increase with increasing A, from values 
lower than ar characteristic for low ice concentrations to values higher than for densely 
packed ice fields (Fig. The relationship between and a suggests the existence of 
different mechanisms dominating cluster formation at low and high ice concentration. In the 
first case, typical clusters consist of the smallest fioes accumulated along the downwind edges 
of the largest ones (Fig. [5^), because of large differences between the respective equilibrium 
velocities. This effect is also present during the initial stages of cluster development at 
higher ice concentrations (not shown), similarly leading to a transient decrease of a. In 
other words, the smallest floes tend to belong to the largest clusters, without contributing 
much to Fc, dominated by the area of the largest floes. Hence, for the largest clusters their 
fgS are hardly larger than the respective radii of their largest members. 
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FIG. 8. Rank-order distributions of effective cluster radius re = (Fc/vr)^/^ in fully clustered states 
for A equal 0.6, 0.725 and 0.9. Black dots show the corresponding distribution of the radii of single 
floes. Values are normalized so that vr^Tg • = 1. 
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FIG. 9. Exponent a of the distributions P{re) of effective cluster radii (equation I16p in function 
of the mean ice concentration A. The horizontal dashed line shows the value of = 1.5. The 
continuous line shows the linear least-square fit to the data for A > 0.7. 

To the contrary, at high ice concentrations, due to high colhsion rates, the fioes tend to 
move with similar translational velocity independently of their size (Fig. [Ht). In this case, 
one large-scale cluster develops, spanning the whole model domain (Figs. and[H]). Only 
the smallest floes can gain enough kinetic energy at collisions so that they are able to escape 
the mega-cluster; and only the smallest fioes find enough empty space in its 'holes' for an 
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FIG. 10. Surface area of the largest cluster Fc,i (relative to the total surface area of all floes, Ftot) 
in function of the mean ice concentration A. For each value of A, the minimum and maximum 
values of Fc^i recorded in the quasi-stationary state is shown. 

independent motion. This results in a > a^. Finding a quantitative explanation for the 
observed dependence a{A), as shown in Fig. [9l is beyond the scope of this paper. 

In terms of the dependence of the largest cluster size on the ice concentration, two regimes 
can be identified, as shown in Fig. [TOl For Aj^O.T, the area of the largest cluster remains 
comparable with the area of the largest floe in the sample (Fig. |8]). For A >~ 0.75, one dom- 
inating mega-cluster develops, occupying over 85% of Ftot- Those two regimes are separated 
by a narrow range of medium, 'critical' ice concentrations Acrit, characterized by a very 
interesting dynamics. Contrary to a relatively smooth development of the largest cluster in 
time, observed for both A ^ Acrit and A ^ Acrit, the transitional range is associated with 
very strong, erratic temporal variations of the size of the largest cluster, corresponding to 
its constant rearrangement, break-ups and re-consolidation. In a sense, when judged by the 
size of the largest cluster, the system never reaches a quasi-stationary state for A ~ Acrit- 

Interestingly, in this transitional state the exponents a characterizing the cluster-size 
distribution P(re) are close to the value of (Figs. [Eland [9]), i.e., the cluster-size distribution 
P(re) mirrors P{r). If a{A) is concerned, the critical range of A separates a region of slowly- 
varying values of a < from a region of fast, approximately linearly increasing values of 
a> ar (Fig. [9]). 

20 



V. DISCUSSION AND CONCLUSIONS 



The model presented in this paper describes a very ideahzed dynamics of MCIZ. Its 
purpose is rather to assist in understanding of basic processes regulating the formation of 
clusters in sea ice, than to reproduce details of any particular real-world situation. Quali- 
tatively, the results of this study reproduce the dominant patterns of motion and clustering 
observed in MCIZ, including realistic granular temperature levels and clusters of small floes 
being pushed in front of a much larger one. Quantitative verification will require, on the 
one hand, experimental data on high-frequency sea-ice motion, unavailable at present, and 
on the other hand, more extensive simulations with a wide range of model parameters {h, 
Ua, e, ar etc.). Possible directions of future research include, but are not limited to: (i) 
formulation of an improved coUisional rheology for MCIZ, based on internal-stress tensor 
components calculated with the SFC model, (ii) better understanding of the atmosphere- 
ice-ocean momentum transfer in MCIZ, resulting from the knowledge of the average motion 
of clustered ice fields in response to winds and currents, (iii) more detailed theoretical anal- 
ysis of phenomena described in this paper, in particular explaining the dependence between 
the exponents of the floe-size and cluster-size distribution, and the nature of the transition 
between the two cluster-size regimes observed at low and high ice concentration, (iv) more 
detailed, in-depth analysis of the analogies between sea ice in MCIZ and other forced poly- 
disperse granular media, and (v) analysis of the role of clustering in atmosphere-ocean heat 
transport and freezing/melting of the ice. 

The last point is particularly interesting in view of a hypothesis formulated recently 
by Toyota and colleagues [itI] who propose that cluster formation (or herding, in their 
nomenclature) may play a role in freezing of neighboring ice fioes, and thus have an influence 
on the fioe-size distribution in MCIZ. Combined with the results of the present study, this 
suggests an interesting possibility of a feedback between FSD and floe clustering. On the 
one hand, the FSD influences processes of cluster formation, including ice concentration 
within clusters and their size distribution. On the other hand, the existence of clusters 
may be important for floe formation in periods with low temperature, contributing to more 
intensive lateral freezing between densely packed fioes. Clusters of fioes frozen together are 
seen in satellite images of MCIZ (see Fig. [11] for an example of a Landsat image of the 
Okhotsk Sea). However, from a single snapshot it is not possible to determine the 'life 
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FIG. 11. Fragment of a Landsat image of the Okhotsk Sea (3. Mar. 2003, 157.1°E, 58.7°N) with 
clusters of thick ice floes consolidated by thin ice. Source: 12l |. 



history' of those floe formations. They could have been produced by freezing within already 
existing clusters, or - alternatively - they could have been a result of a divergence of an 
initially densely packed, frozen floe fleld and of subsequent random breaking of thin ice 
occupying spaces between thick floes. In the second case, the end result could have formed 
without the clustering-freezing-FSD feedback proposed above. Estimating the plausibility 
of that feedback is not possible without more observational data. 

Even though substantial progress has been made in recent years in observational tech- 
niques of the polar environment, including sea ice, MCIZ is still a very demanding, hard-to- 
explore medium. Consequently, even though it plays a crucial role in seasonal and long-term 
expansion and retreat of the ice cover in polar seas - one of important proxies of the climate 
change - the dynamics of MCIZ still remains poorly understood. Numerical models, like 
the one developed in this study, provide a valuable tool enabling to improve our knowledge 
and to gain more insight into functioning of that complex and fascinating environment. 



Appendix A: Rank-ordering statistics and maximum-likelihood estimation 



To avoid drawbacks of binning in the analysis of distributions of the modeled quantities, 
rank-ordering technique is used throughout this paper 
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. It is based on data sorted in 



descending order and - graphically - plotted versus the ranks (as in Fig. |8]). 
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1. Calculation of the floe diameters 



It is assumed that the floe-size probabihty density function (pdf ) P{d) is given by a power 
law P{d) = Crf^^^"'' with d > and a,. > 1, so that the mean d exists. For prescribed a^, 
d and the flnite sample size A^, acting as a constraint, the maximum-likelihood estimation 
gives the most probable value of the nth of N elements as [371, p. 171]: 

l/ar 



arN + 1 
an + 1 



C 



(Al) 



The influence of the flnite sample size is visible in the range of the largest elements in the 
sample: with decreasing A^, the probability of encountering very large values decreases. This 
effect is seen in Fig. |8] as a deviation of the rank-order distribution of floe diameters from 
a straight line. 

In the SFC model, formula (lAll) is used with C = 1, and the obtained values are subse- 
quently re-scaled so that X]^=i = Nd. For a given ice concentration A, the length L of 
the computational domain is calculated as L = (7r/(4yl) Ylin=i ^n)^^'^- 



2. Estimation of power-law exponents 

Let P{x) denote a power-law probability density function of variable X, represented by 
A^ rank-ordered data points x„ {n = 1, . . . , N). A maximum-likelihood estimation of the 
exponent a of P{x), based on the sub-sample (xkjX^), where (1 < A; < A^), is given by js^]: 

(A2) 

The standard deviation Aa^ of this estimate decreases with increasing sample size: Aa^ = 
ak{N - k + 



ak = N 



N 



n=k 



xn 
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